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1. Introduction 



Consider the simple communication system model 
shown in Fig. 1. The goal is to transmit one of M 
possible symbols, i.e., an M-ary signaling system, over 
a memoryless additive noise channel. We will assume 
all signals are discrete-time with T samples. The trans- 
mitter assigns a unique signal ^;„ : { 1 , . . . , T } — > [R to 
each symbol m G { 1 , . . . , M } . It is this signal that is 
sent through the channel. At the other end, the received 
signal is 

y\t^=Sm\t^ + n\t\ ?=l,...,r, 

where f2:{l,...,r}— >[Risa noise process, and the 
job of the receiver is to decide which symbol was trans- 
mitted. Our goal is to design a set of signals Sm, 
m = 1, . . . , M, which maximize, subject to constraints 
on the signals, the probability of a correct decision by 
the receiver given a particular channel noise distribu- 
tion. 



Of course, in order to design an optimal signal set, the 
action of the channel and the receiver must be com- 
pletely specified. For the channel, we assume the noise 
process is independent and identically distributed (iid) 
with distribution ;7iv. Further, we assume that the noise 
process is independent of the symbol being transmitted. 
Our detection problem falls into the class of M-ary 
Bayesian hypothesis testing problems where, for 
m = 1, . . . , M, the hypotheses are defined as follows, 

Hm : y\t^ = s^[t] -\-n[tl ? = 1, . . . , r. 

To simplify notation, define the received signal vector 

Finally, it is assumed that the receiver was designed 
using the minimum average probability of error crite- 
rion (or the uniform cost criterion). It is well known that 
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(see, e.g., Sec. IV.B of Ref. [20]), under our assump- 
tions, the optimal receiver is the maximum a posteriori 
probability (MAP) detector. Specifically, the optimal 
receiver chooses 



m(y) = arg m20i{p(HJy) I m = 1, . 



,M} 



i.e., the hypothesis with the largest probability given the 
observation y. 

Clearly, the receiver will make an error if hypothesis 
Hm is true, but 

p(H^-\y)>p(HJyl 

for some m' ^ m. Thus, the probability of a correct 
decision under hypothesis Hfn is 

/7({correct decision} I H^) = p([p(HJy) > p(Hnr-\y), 

\fm'i^m} \HJ 



-({■%i^><'- '■"■*''] 



H, 



where, in order to put things in terms of the familiar 
log-likelihood ratio, we have assumed p(H^'\y) > for 
all y, m' G { 1, . . . , M}. For the signal set design prob- 
lem considered here, no knowledge of the prior distribu- 
tion on the hypotheses //^, m = 1, . . . , M, will be as- 
sumed. Of course, the conditional distribution /? (//^ I y) 
is known since, given a signal set, this distribution is 
completely determined by the distribution on the chan- 
nel noise. Specifically, in view of our assumptions, 

T 

p(HJy) = ^PN(y[t] -s^[t]), 

t=i 

If the prior distribution were known, the quantity to 
be maximized could be expanded as 



As p{Hm) is not assumed to be known, the worst-case 
prior distribution will be used to compute /?({ correct 
decision}) for any given signal set. In particular, let 



1, y^^O, m=l,...,M 



The goal will be to find signal sets which maximize 

min A /?({ correct decision} I H^) * 7m. 

It is not difficult to show that this is equivalent to max- 
imizing 



min /?({ correct decision} I //;„). 

me{l,...,M} 



(1) 



A standard assumption in transmitter design is that 
the signals are restricted to be of the form 



^mM = E«"''^^W. 



(2) 



where c^^ : { 1 , . . . , T} ^ [R, /: = 1 , . . . , ^, are given 
basis functions and a""'^ G [R, m = 1, . . . , M, 
/:= 1, . . . , /T, are the free parameters. Finally, due to 
power limitations in the transmitter, the signals are 
forced to satisfy some type of power constraint, either 
peak amplitude or average energy. In this paper, we will 
assume a peak amplitude constraint, i.e., 

I s^[t]\ < C, m = 1, . . . , M, r = 1 , . . . , r, (3) 

where C > is given. Note that we could just as easily 
have considered an average energy constraint in our 
formulation. Our design problem is thus reduced to 
choosing parameters a""'^ in order to maximize Eq. (1), 
subject to the constraints Eq. (3). 



/?({correctdecision})=2/? ({correct decision} \Hj,,)'p{Hm\ 
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2. The Optimization Problem 

The details of casting the design problem introduced 
in the previous section in such a way that it may be 
solved using standard nonlinear programming al- 
gorithms are presented in this section. The design of 
optimal signal sets under the assumption of Gaussian 
noise has been well studied (see, e.g., Ref. [22]). In fact, 
a gradient-based first-order algorithm was developed 
and analyzed in Ref. [5] for the case of Gaussian noise, 
K=2 basis functions, and an average energy constraint 
on the signals. The performance of optimal detectors in 
the presence of non-Gaussian noise as a function of the 
signal set was first studied by Johnson and Orsak in Ref. 
[10]. It was shown in Ref. [10] that the dependence of 
detector performance on the signal set is related to the 
Kullback-Leibler (KL) distance between distributions 
for the various hypotheses. Based on this work, Gocken- 
bach and Kearsley Ref. [7] proposed the nonlinear pro- 
gramming (NLP) formulation of the signal set design 
problem which is considered here. 

Given our assumptions on the noise process, the log- 
likelihood ratio may be written 



In 



pjHJy) _y p{HJy\i]) 

PiHJly) tt PiH^'\y[t])' 



Note that, since randomness only enters the received 
signal through the additive noise process, when hypoth- 
esis Htn is true, the receiver computes 

p(HJy[t])=p^(n[t]), 

and, for m' i^ m, 

p(Hm'\y[t])=PN(n[t] + (s^'[t] - s^[t])). 

Thus, upon substitution, the statistic of interest to us is 

1 PJHJy) ^ y . ( PNJnit]) 

p(HJy) tr W(^M + (^m'[n -^.[?]))/* 

Now, assuming the variance of the statistic (4) does 
not change as we vary m' i^ m, maximizing ;?({ correct 
decision} \Hfn)is equivalent to maximizing the expected 
value of the statistic Eq. (4) for each m' i^ m. That is, 
under hypothesis Hm , the probabihty of correctly choos- 
ing Hm is maximized if we maximize 



min^lilnf , ,,/f\'}^ 



sAm 



H„ 



Thus, for the signal design problem considered here, we 
may equivalently use 



mm mm 

mG{l M} m'^m L ;^i 



It! \PN(n[t] + (s„.[t]-sAmJ 



Hm 



as the objective function. Define the function 



i.e., the KL distance between the noise distribution and 
the noise distribution shifted by —5. Note that if we 
assume a symmetric distribution for the noise (this is not 
a restrictive assumption), then Kj^{') will be an even 
function. It is not difficult to show that 



M) 

vt^-sAm 



[ft \PN{n[t] -\- {Sm' 
T 



Hm 



Define 
a = {a 



. , a 



. , a 



M,\ 



, , a^'^) - 



Substituting the expansion Eq. (2), we see that, under 
our assumptions, the signal set design problem is equiv- 
alent to solving the optimization problem 



mm max 



^ t=\ \k=\ / 



m,m 



:{1, . . . ,M}, m'>m 



s.t. (2] «"''(/>/: [^]) ^C, m=l,...,M, r=l,...,r. 

(SS) 

It is only necessary to consider m' > m since ^iv(*) is an 
even function. 

The problem (SS) is already in a form which may be 
handled by algorithms which tackle inequahty con- 
strained mini-max problems. Such algorithms have been 
developed by, e.g., Kiwiel in Ref. [13], Panier and Tits 
in Ref. [16], and Zhou in Ref. [25], all in the context of 
feasible iterates. In Ref. [25], Zhou extends the non- 
monotone line search-based algorithm of Ref. [26] to 
handle the constrained case. The algorithm of Ref. [16] 
extends the feasible SQP algorithm of Ref. [17] to han- 
dle mini-max objective functions. A recent algorithm for 
the constrained mini-max problem which does not gen- 
erate feasible iterates is the augmented Lagrangian ap- 
proach of Rustem and Nguyen Ref. [23]. The problem 
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may be converted into an equivalent single objective 
smooth nonlinear programming problem by adding an 
additional variable. This is the approach taken by Gock- 
enbach and Kearsley in Ref. [7]. Additionally, they add 
a "regularization" term (possibly zero) to the converted 
objective. Specifically, consider 



mm 



y'- €r\\a\\l 



s.t. 



t=\ ^k=\ 

y^ m', m G {1, . . . ,M}, 



2«"'^c^fe[^] ^C", m= 1,...,M, r=l,...,r, 

y^o, 

(^^') 



where e^ is small (possibly zero). 

It turns out that problems (SS^ and {SS'^ are difficult 
to solve using standard off-the-shelf nonlinear program- 
ming codes. To begin with, it was observed in Ref. [7] 
that outside of the feasible region, the linearized con- 
straints for problem {SS) are often inconsistent, i.e. no 
feasible solution exists. This can be a big problem for 
sequential quadratic programming (SQP) based al- 
gorithms, generally considered the most successful al- 
gorithms available for NLP problems with a reasonable 
number of variables. Of course, with feasible iterates, 
the hnearized constraints are always consistent and the 
solutions of the QP sub-problems are always well-de- 
fined. As an alternative to feasible iterates, there are 
infeasible SQP-based algorithms which use special 
techniques to deal with inconsistent QPs (see, e.g., Ref. 
[24, 12, 4]). Another difficulty in applying a local NLP 
algorithm is that problem {SS) has many local solutions 
which may prevent convergence to a global solution. 



3. Local Algorithms 

Sequential Quadratic Programming (SQP) has 
evolved into a broad classification encompassing a vari- 
ety of algorithms. When the number of variables is not 
too large, SQP algorithms are widely acknowledged to 
be the most successful algorithms available for solving 
smooth nonhnear programming problems. For an excel- 
lent recent survey of SQP algorithms, and the theory 
behind them, see Ref. [4]. 



In general, an SQP algorithm is characterized as one 
in which a quadratic model of the problem is formed at 
the current estimate of the solution and is solved in order 
to construct the next estimate of the solution. Typically, 
in order to ensure global convergence, a suitable merit 
function is used to perform a line search in the direction 
provided by the solution of the quadratic model. While 
such algorithms are potentially very fast, the local rate 
of convergence is critically dependent upon the type of 
second order information utilized in the quadratic model 
as well as the method by which this information is 
updated. 

3.1 Infeasible SQP 

Numerous SQP algorithms that do not require iterates 
to remain feasible have been suggested by researchers 
(e.g., Refs. [3, 6, 24] among others). Because of the 
nonlinear nature of the constaints appearing in this 
specific class of problems, the linearizations employed 
by SQP are frequently inconsistent. As a result, con- 
straint perturbation techniques must be employed to en- 
sure that the algorithm can generate a descent direction 
at every iteration. This, at least in part, explains the 
superior performance of hybrid SQP algorithms re- 
ported on here (see Tables 2 and 3). These algorithms 
are particularly desirable because they do not require a 
feasible starting point. 

3.2 Feasible SQP 

In Ref. [19, 17, 18, 9], variations on the standard SQP 
iteration are proposed which generate iterates lying 
within the feasible set. Such methods are sometimes 
referred to as "Feasible SQP" (or FSQP) algorithms. It 
was observed that requiring feasible iterates has both 
algorithmic and apphcation-oriented advantages. Al- 
gorithmically, feasible iterates are desirable because 

• The Quadratic Programming (QP) subproblems are 
always consistent, i.e., a feasible solution always ex- 
ists, and 

• The objective function may be used directly as a merit 
function in the hne search. 

State of the art SQP algorithms typically include com- 
plex schemes to deal with inconsistent QPs. Further, the 
choice of an appropriate merit function (to enforce 
global convergence) is not always clear. Requiring feasi- 
ble iterates ehminates these issues. In order to overcome 
the problem of inconsistent QPs in this work, we use the 
CFSQP implementation Ref. [14] of the FSQP al- 
gorithm proposed and analyzed in Ref. [18] as our local 
algorithm. 
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4. Global Solutions 



4.1 Stochastic Approach 

In one attempt to overcome the problem of local solu- 
tions, we use a stochastic two-phase method (see, e.g., 
Ref. [1]) where random initial points are generated in 
the feasible region and a local method is repeatedly 
applied to a subset of these points. Such an approach 
may be thought of as simply a "smart" way of generat- 
ing many initial points for our fast local algorithm with 
the hopes of eventually identifying a global solution. 
Specifically, we will use the Multi-Level Single Linkage 
(MLSL) approach Ref. [1, 11], which is known to find, 
with probability one, all local minima (hence the global 
minima) in a finite number of iterations. 

Let M denote the cumulative set of local minimizers 
identified by the MLSL algorithm. At iteration €, for 
some integer A'^ > fixed, we generate A'^ points 
Q;(/-i)Ar+i, . . . , a^N distributed uniformly over the 
feasible set X for (SS), For each of the points 
«/ G {«],..., a^} we check to see if there is another 
point aj within a "critical distance" r^ of «/ which also 
has a smaller objective value. If not, then the local al- 
gorithm, call it i£ is applied with initial point «/ and the 
computed local minimizer is added to the set M, After 
all points are checked, r^is updated, ^is set to ^+ 1 and 
the process is repeated. At any given iteration, the local 
maximizer with the smallest objective value is our cur- 
rent estimate of the global solution. For ease of notation, 
define the (mini-max) objective 



r T I K 

F{a) ^ max -^ kJTZ («'"''' " «"^'')</>^[n 



m,m 



{1, . . . ,M}, m'>m 



Further, let ££(«) denote the local minimizer obtained 
when a local algorithm is applied to problem (SS) with 
initial point a. The following algorithm statement is 
adapted from Ref. [1]. 



Step J . set / ^ / -I- 1 . 

if i < €N then goto Step 2 . 

else set € <^ €-\- I and goto Step 1 . 

It remains to specify how we select the critical dis- 
tance r^, the definition of the metric 11*11^ we use for 
signal sets (as parameterized by a ), and how we gener- 
ate the sample points. Following Ref. [1], we use 



^/ = 



1 



(r(i 



-F 71/2) ■ 



,<x,.i^)' 



where n is the number of variables (MK for our prob- 
lem), m(X) is the volume of the feasible region, and 
^ > 2. To compute m(X), note that in view of symmetry 
with respect to the signals, 

m(X)=A, 

where A is the volume of the feasible region for the 
parameters corresponding to one signal (recall, M is the 
number of signals). The quantity A is easily estimated 
using a Monte Carlo technique. 

Note that, for our problem, as far as the transmitter is 
concerned, a given signal set is unchanged if we were to 
swap the coefficients q;'"^'^, A:= 1, . . . ,^, with a"'^'\ 
/:= 1, . . . , ^, where rrii + m^. The distance "metric" 
we use in Algorithm MLSL should take this symmetry 
into account. Consider the following procedure for com- 
puting the distance between signal sets parameterized 
by ax and a^, 

function dist(ai, 0:2) ( 

for /=1,...,M- 1 do { 



a^if 



f; = min j^ {a\ 



ji = arg min 



yE{l,...,M}/{ji,...,j--i} 



a^'f 



jE{l,...,M}/{j 



/i, ...,jVi}| 



return 



{^4 



I 



Algorithm MLSL 

Step 0, set € ^l,M^0, 

Step 7. generate A^ points Q;(^-i)Ar+i, 
formly over X. 
set f ^ 1. 



, am urn- 



Step 2 . if (3j s.t. F(aj) < F{ai) and Wat - ajM < re) 
then goto Step 3 . 

else set il ^ il U {^{ad}. 



This is not a metric in the strict sense of the definition, 
though it suffices for our purposes and we will set 

\\ci\ — ciiWs = dist(Q;i,Q;2). 

To aid the generation of sample points, before starting 
the MLSL loop we compute the smallest box which 
contains the feasible setX. By symmetry with respect to 
the signals, we can do this by solving 2K hnear pro- 
grams. Specifically, let a^ G [R, /: = 1 , . . . , ^ be defined 
as the optimal value of the hnear program (LP) 
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max a 



s.t. 2«'''<Afe[^] ^c^ r=i,...,r, 



9=1 
K 



(U,) 



q=\ 



and let a^ E [R, ^ = 1, . . . , ^ be defined as the optimal 
value of the LP 



Thus, the problem has grown from KM to KM(M — 1) 
variables. Let a denote the vector of all such coeffi- 
cients. As mentioned above, the coefficients a^'^'^ and 
a^'"'^, ^ = 1 , . . . , /T will be used only to compute the 
interaction between signals m and €, Correspondingly, 
define the objective function which encapsulates this 
interaction as 

tei \t=i / 



mm a 



a ^',...,a ^ 



S.t. ^a'^'(l)k[t] ^c, r=i,...,r. 



9=1 
K 



iU) 



q=\ 

Then, it should be clear that 

Z C ^ = {a G [R^^la"^'^ G [a^ a:^], 

^=l,...,^,m=l,...,M}. 

Using standard random number generators, it is a simple 
matter to choose samples from the uniform distribution 
on the box S3. Thus, for Step 1 of Algorithm MLSL, we 
repeatedly generate samples from the uniform distribu- 
tion on S3, discarding those which do not lie in Z, until 
we find A^ which do lie in Z. It should be clear that such 
a procedure is equivalent to drawing A^ samples from the 
uniform distribution on Z. 

4.2 Expanded Space Approach 

In this section we describe a reformulation of the 
problem along the lines of that considered in Ref. [8] in 
the context of molecular conformation. The essence of 
the approach is to add variables in such a way that the 
global solution of the new equivalent problem has an 
enlarged basin of attraction. While the method does not 
guarantee convergence to a global solution, it does in- 
crease the likelihood. 

To use such an approach in our context, we assign to 
each signal M — \ sets of "auxiliary" coefficients. Each 
set will be used exclusively for computing the interac- 
tion with one particular other signal. For a given signal, 
the method will asymptotically force these sets of auxil- 
iary coefficients to coalesce, as they must in order for 
the solution to be feasible for the original problem. Let 
the auxiliary coefficients be denoted a""'^'^, where m G 
{l,...,M},yfG{l,...,M}/{m},and^G{l,...,^}. 



A schematic representation of the auxiliary coefficients 
and their interactions is given in Fig. 2 for the case 
M=3. 

The constraint that the sets of auxiliary coefficients 
for each signal must all be equal at a feasible solution is 
easily expressed as a set of linear equality constraints. In 
particular, it is not difficult to show that there exists a 
linear operator W, a KM{M - 2) X KM{M - 1) ma- 
trix, whose kernel is precisely the set of all vectors 
satisfying this constraint. Finally, in view of this con- 
straint, it is only necessary to enforce the power con- 
straint on one set of auxiliary coefficients for each sig- 
nal. Thus, our "equivalent" expanded space problem is 

min maxlK^^'fa) I m, m' G { I, . . . , Ml, m' > m 1 

s.t. ('2a'-^^4>dt]] ^ C\ t=\,...,T, 

(ES) 

i^dr^'^'Miyj ^c\ r = I, . . . , r, m = 2, . . . , M, 

^fc=l / 

Following Ref. [8], instead of attempting to solve this 
problem directly, we incorporate the equality constraints 
into the objective functions as a quadratic penalty term. 
This allows us to approach solutions from "infeasible" 
points by carefully increasing the penalty parameter and 
asymptotically forcing the auxiliary coefficients to coa- 
lesce for each signal. Specifically, define the penalized 
objective 



P..,Aa; p) ^ K^,^-(a) + i p^a^W^Wa, 



for m,m' G {1, . . . ,M}, m'> m. For a fixed value of 
p > 0, we can solve (using, e.g., a mini-max SQP-type 
algorithm) the problem 
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min max{P;„,;„'(a;p) I m,m'E{l,..., M},m'>m} 



-1) 

s.t. (i«'•'•*(^.[f]ysc^ t=\,...,T, 

PESip) 
2a'"•'■V,[f]y<C^ t=\ T, m = 2 M. 

k=\ / 

Let S(a°, p) denote the solution of PES(p) obtained 
using a local algorithm (such as those discussed in Sec. 
3) starting with the initial point oP. Using the solution 
just obtained as the next initial point, the process may be 
repeated after appropriately increasing the penalty 
parameter p. At any time, a candidate global solution for 
the original problem (SS) may be computed by 
"projecting" a solution a = Z{a\ p) of PES{p) into W"" 
and solving (SS) using a local algorithm with the pro- 
jected vector as the initial point. We will denote the 
projection operation Proj(a) and, using the notation 
from Sec. 4.1, the computed local solution for (SS) 
starting from the initial point Proj(a) will be denoted 
S(Proj(a)). One possible method to compute the re- 
quired projection is to simply take the arithmetic average 
of the corresponding components of the auxiliary coef- 
ficients, i.e., if a = Proj(a), then componentwise 



M - 1 £l "^ ' 



m=l,...,M,)t= 1,...,^. 

It remains to specify how we update the penalty 
parameter p. At "major" iteration^ /, the penalty 
parameter will be increased by a multiplicative factor, 
call it 8i, i.e., 

p,+i = 8i ' Pi. 

In addition, we will decrease the factor 5, when the 
projection for the current major iteration produces an 
estimate which does not improve upon the previous esti- 
mate. If Pi is increased too fast for a given problem we 
could get trapped in a local solution. The precise al- 
gorithm statement follows. 

Algorithm ES 

Data, oo e [R^^t^-i), p, > 0, 5o > 1. 
Parameters . €>0. 



A major iteration is defined here as on solution of PES(p) and one 
update of p . 
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Step 0, set / <— 1, ao <— Proj(5;o). 

Step 1 . compute a, = S(a/-i, p/). 

Step 2. compute a, = S(Proj(a/)). 

Step J. if (F(a,) > F(a,-i) + 6) then set 

else set 6^ <— S^-i. 

5?^;? 4. set P/+1 <— 6/ * p,, f <— / + 1. 

5?^;? 5 . goto 5?^;? 7 . 

4.3 Homotopy Approach 

Suppose that for a given set of basis functions (and 
fixed problem size) we know the globally optimal signal 
set for one particular noise distribution, say p\{'). For 
example, in many cases it is possible to analytically 
compute the globally optimal signal set for Gaussian 
noise. In this section, we discuss a method for comput- 
ing candidate globally optimal signal sets for a different 
noise distribution, say pii'), based on this knowledge. 

The so-called homotopy approach relies upon the 
observation that, for AG [0,1 ], 

;?n(t; A) = (1 - X)pi{t) -h Xpiir) 

is a valid noise distribution. Gradually increasing A and 
repeatedly applying a local algorithm, the computed 
optimal signal set should trace out a continuous "path" 
from that for px to that for p2> Let a(A , ao ) denote the 
computed optimal solution of {SS) obtained via a local 
algorithm starting with the initial point ao and using the 
noise distribution p^ (* ; A). At iteration f, we compute 

a/+i = a{Xi, a,), 

and the homotopy parameter A^ is updated. The proper 
updating of A/ is critical. Clearly, we should have 

hm Ai = 1 . 

Further, the rate of convergence of the method is di- 
rectly related to how quickly A/ is increased. On the 
other hand, if A/ is increased too quickly then a/+i may 
"jump" to a new path of minimizers. For prehminary 
tests we simply increment A/ by a small, fixed, amount. 

Algorithm HOM 

Data, aoG [R''^0<6a«1. 
Step 0. set i ^ 1, Ao ^ 6a. 



Step 7. compute a, <— a(A/_i, a/-i). 

Step 2 . if A/ = 1 then stop. 

Step J. set A/+1 <— min{A/-i- 6a, 1}. 

Step 4 . set / <— / -H 1 and goto Step 1 . 

Numerically, one of the biggest challenges associated 
with this approach is the computation of the KL distance 
K^ and the derivative of the KL distance K'^, For the 
distributions considered in this work these functions are 
readily obtained analytically. Unfortunately, this is no 
longer possible for convex combinations of these distri- 
butions as considered in this section. Consequently, we 
are forced to turn to numerical quadrature. For the pre- 
liminary numerical implementations we use a simple 
infinite trapezoid rule, i.e. we use the approximation 

For those integrands with a slowly decaying tail we use 
the change of variables 



T(0 = 



5. Numerical Results 

Following Ref. [7], we consider the noise distribu- 
tions Pn listed in Table 1 . For the definition of the Gen- 
eralized Gaussian distribution, the constant a is defined 

as 



a^ 



7T(i/4) y 
r(3/4) I 



For our numerical experiments, we assume cr= 1. The 
case ^ = 2 is of common interest, and we use either a 
sine-sine basis 

I y| sin(27TWir/r), y| sin(27TW2?/r)|, 

or a sine-cosine basis 

I y- sin(27TWir/r), y- cos{2iio},tlT) I, 

where W] = 10 and 0*2 = 1 L When K=2wq can display 
the results in the plane as familiar signal constellations . 
Finally, we run experiments for M=8, 16 signals, 
T=50 samples, and with an amplitude bound of 
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Table 1. Noise distributions and the associated KL distance function 



Noise distribution 



PN(r) 



Kn(8) 



Gaussian 
Laplacian 
Hyperbolic Secant 
Generalized Gaussian 
Cauchy 



^p(^) 



I 






<=^) 



— sech -— ] 
a \2(t; 



1 



2T(5/4)a 
1 



exp 



a } 



TTO-Cl + Cr/o-)") 



_5_ 
la" 

V2-I6I /-V2-I6l\ , 
+ expl I - 1 



/-V2 ■ I6I\ 



r^(3/4) 
r'(i/4) 

'"(■*i>) 






C = VIO. Note that, for M= 16, problem (5'5') has 32 
variables, 120 objective functions, and 800 constraints. 

In order to demonstrate the need for special-purpose 
SQP codes, we attempted to solve the problem using 
VF02AD^ from the Harwell subroutine library Ref. 
[15], a standard SQP code based on Powell's algorithm 
Ref. [21] and a hybrid SQP code recently developed Ref. 
[3] and analyzed Ref. [2]. These codes do not directly 
solve mini-max problems, we used the formulation (SS'^ 
suggested in Ref. [7]. In Tables 2 and 3, we list the 
number of times VF02AD and the BKT SQP algorithm 
Ref. [3] successfully converged to a local minimizer out 
of 20 trials for a given noise distribution and basis (and 
regularization parameter). In parenthesis we report on 
the number of times each algorithm succeeded in con- 
verging to the global minimizer. For each trial the initial 
point was drawn from the uniform distribution over the 
feasible set. 

It is clear from the table that the standard SQP al- 
gorithm had little success converging to a local solution. 
The failures were essentially always due to inconsistent 
constraints in the QP sub-problem. In our trials, neither 
the Feasible SQP algorithm nor the hybrid BKT SQP 
algorithm failed to converge to, at least, a local solution. 

We ran Algorithm MLSL for 20 different problem 
instances. The algorithm was stopped after it appeared 
that no better local minimizers would be found (i.e., the 
estimate of the global minimum remained constant for 



Table 2. Number of successful runs for VF02AD out of 20 trials 



^ Certain commercial equipment, instruments, or materials are identi- 
fied in this paper to foster understanding. Such identification does not 
imply recommendation or endorsement by the National Institute of 
Standards and Technology, nor does it imply that the materials or 
equipment identified are necessarily the best available for the purpose. 



Noise distribution 


sine-sine 


sine-cosine 


sine-cosine 




(€, = 0) 


(€.= 


0) 


(€.= 10-^) 


Gaussian 


4(1) 







1(0) 


Laplacian 


6(1) 







1(0) 


Hyperbolic Secant 


5(1) 










Generalized Gaussian 


6(2) 










Cauchy 


2(1) 











Table 3. Number of successful runs of BKT SQP algorithm out of 20 
trials 



Noise distribution 


sine-sine 


sine-cosine 


sine-cosine 




(€, = 0) 


(€,= 0) 


(€.= 10-^) 


Gaussian 


20(8) 


18(2) 


18(1) 


Laplacian 


20(4) 


20(3) 


19(5) 


Hyperbolic Secant 


20(5) 


20(1) 


20(4) 


Generalized Gaussian 


20(4) 


18(1) 


20(1) 


Cauchy 


20(2) 


19(1) 


20(3) 



several MLSL iterations). In Tables 4 and 5 we hst our 
computed minimum values for instances of {SS) with 
M = 8 and M = 16, respectively. Note that our solutions 
respectively. Note that our solutions agree with those 
reported in Ref. [7]. In all cases, execution was termi- 
nated after no more than 10 to 15 minutes. In Figs. 3 
through 8 we show the optimal signal constellations for 
several of the instances of {SS) corresponding to the 
optimal values listed in Tables 4 and 5. 



449 



Volume 106, Number 2, March-April 2001 

Journal of Research of the National Institute of Standards and Technology 



Table 4. Optimal computed values for signal set design with M =S 



Table 5. Optimal computed values for signal set design with M = 16 



Noise 


Basis 


Fia*) 


Noise 


Basis 


Fia*) 


Gaussian 


sine-sine 


-69.793 


Gaussian 


sine- sine 


-29.314 




sine-cosine 


-97.551 




sine-cosine 


-39.742 


Laplacian 


sine-sine 


-63.122 


Laplacian 


sine- sine 


-32.370 




sine-cosine 


-84.463 




sine-cosine 


-44.076 


Hyperbolic Secant 


sine-sine 


-61.093 


Hyperbolic Secant 


sine -sine 


-29.577 




sine-cosine 


-83.196 




sine-cosine 


-40.500 


Generalized Gaussian 


sine -sine 


-189.09 


Generalized Gaussian 


sine -sine 


-57.829 




sine-cosine 


-264.18 




sine-cosine 


-76.138 


Cauchy 


sine-sine 


-22.731 


Cauchy 


sine -sine 


-11.426 




sine-cosine 


-30.673 




sine-cosine 


-15.688 




6. Conclusions 

In this paper we have presented an important and 
difficult optimization problem along with a broad arse- 
nal of numerical optimization algorithms and modern 
enhancements that can be used to solve it. These prob- 
lems are not "large-scale" by modem computing stan- 
dards but they are, nonetheless, extremely difficult 
problems to solve efficiently. 

Numerical solutions to these problems were located 
using SQP methods embedded into stochastic global 



algorithms. Numerous numerical tests suggest that this 
embedding procedure can significantly improve the 
performance of the SQP method. 

Because there are so many different algorithms and 
implementations for the solution of nonlinear program- 
ming problems there is a need to create an accepted 
collection of test problems arising from the application 
of optimization. Because of the difficulties it poses, this 
family of problems is a natural candidate to use as a 
practical and significant test problem. 
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